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INTRODUCTION 


The  objective  of  this  study  contract  is  to  develop  a  practi- 

r  ^ 

cal  single  mode  bending  loss  model  for  special  .  fibers  critical  to 
several  future  Army  Weapon  Systems.  The  model  will  facilitate  the 
selection  of  fiber  and  aid  the  design  of  high  speed  missile  payout 
canisters  used  in  major  Army  fiber  optics  systems  such  as  FOG-M 
and  AAWS-M.  The  initial  effort  will  be  directed  to  study  various 
bending  induced  loss  mechanisms  in  fiber.  A  theoretical  bending 
loss  model,  expressed  in  appropriate  computer  algorithms  is  being 
formulated.  Practical  fiber  characterization  schemes  will  be  de¬ 
vised  to  yield  relevant  input  data  to  the  loss  model.  The  model 
will  then  be  modified  to  improve  its  adequacy  for  bending  loss 
analysis.  Environmental  effects  on  fiber  bending  loss  will  be  in¬ 
vestigated.  Reduction  of  temperature  induced  fiber  loss  of  mis¬ 
sile  payout  bobbins  and  field  deployable  fiber  cables  is  the  ulti¬ 
mate  goal . 


PROGRESS 

In  the  first  phase  of  the  Basic  Program  we  have  laid  the 
foundation  for  the  real  thrust  of  the  project.  The  schedule  is 
shown  in  Figure  1.  An  oral  progress  report  was  given  to  ARO  and 
CECOM  personnel  at  CECOM  on  September  26,  1986  and  represents  the 
detailed  portion  of  this  interim  report.  Progress  is  summarized 
and  documented  in  this  report.  A  copy  of  the  oral  presentation  is 
shown  as  Appendix  A.  A  literature  survey  was  conducted 
and  completed  on  the  subject  of  single  mode  fiber  theory  and 
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Figure  1.  Basic  Program  Schedule 


bending  loss  phenomena.  The  survey  has  produced  a  reference  list 
of  76  titles,  presented  as  Appendix  B.  It  shows  extensive  work  to 
date  in  pursuit  of  fundamental  understanding  of  fiber  bending 


loss.  Numerous  articles  have  been  published  on  both  the  macrobend 


and  microbend  fiber  loss  study.  However,  many  experimental  data 
still  can  not  be  fully  explained  by  the  existing  theory.  There  is 
no  unified  theoretical  equation  or  a  single  model  that  adequately 
predict  actual  fiber  bending  loss.  In  addition,  there  are  prob¬ 
lems  generated  by  the  different  analytical  approaches  employed 
during  fiber  bending  loss  research.  Our  immediate  effort  is  to 
review  these  existing  theories  and  to  formulate  a  comprehensive 
single  mode  fiber  loss  model  that  combines  the  output  of  past  re- 
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search  efforts  with  new  work.  The  fiber  parameters  and  measured 
bend  loss  are  the  inputs  and  the  output  of  the  mode. 


The  basic  mechanism  of  bend  induced  loss  on  a  single  mode  fi¬ 
ber  is  a  mode  coupling  process  taking  place  between  a  guided  mode 
(HE11)  and  the  radiation  modes  of  fiber.  Specifically,  the  ra¬ 
diation  modes  include  both  cladding  modes  and  air  modes.  A  mode 
coupling  into  the  cladding  modes,  in  which  the  optical  power  is 
still  trapped  in  the  fiber  cladding,  often  creates  a  slow  power 
leakage  along  the  fiber  length.  A  mode  coupling  to  the  air  modes 
will  cause  a  radiation  loss  as  the  optical  power  actually  radiates 
out  the  fiber.  The  fiber  bend  loss  mechanism  can  be  roughly 
divided  into  two  categories,  depending  on  its  physical  dimensions 
and  the  abruptness  of  the  bend.  The  categories  are  shown  in  Fig¬ 
ure  2 . 

•  MACROBENDS:  LARGE  BENDING  CURVATUKt  ( R  >>  b)  THAT  COVERS  A  FIBER  LENGTH  OF  L  >;>  b 

SPOOL;  MANDREL 
CORRUGATION  PLATE 
PAYOUT  PEEL  POINT 

l 

•  MICROBENDS:  SHARP  BENDING  CURVATURE  <R~  b)  THAT  INVOLVES  LOCAL  AXIAL  DISPLACEMENTS 

OF  A  FEW  MICRONMETERS  AND  SPATIAL  WAVELENGTHS  OF  A  FEW  MILLIMETERS. 

FIBER  DIMENSION  NONUNIFORMITY  (CORE  SIZE) 

CABLING  STRESS,  PACKING  PRESSURE 
BOBBIN  WINDING  STRESS 

_  ENVIRONMENTAL  STRESS  UN  BOBBIN  WOUND  FIBER  (LOW  TEMPERATURE) 


Macrobend  generally  refers  to  a  bend  curvature  several  orders 


larger  than  the  optical  wavelength.  It  can  introduce  radiation 
loss  as  the  result  of  field  deformation  on  the  fiber  guided  mode. 
The  radiation  loss  depends  strongly  on  the  bend  curvature  and  the 
fiber  index  profile.  In  contrast,  microbend  refers  to  the  micro¬ 
scopic  random  deviation  of  fiber  axis  from  its  natural  straight 
condition  defined  by  the  original  drawing  of  the  fiber.  It  can  be 
introduced  on  fiber  by  cabling,  winding,  and  ambient  environment 
change.  Microbends  often  generate  gentle  mode  coupling  between 
the  guided  mode  and  the  cladding  modes  of  fiber  and  generally  lead 
to  a  small  optical  power  loss  in  fiber  over  a  long  length.  The 
microbending  loss  is  known  to  depend  on  fiber  structure,  jacketing 
material,  cabling  design,  winding  condition,  and  ambient  condi¬ 
tions.  In  theory,  microbending  loss  is  a  complex  process  that  of¬ 
ten  requires  statistical  methodology  to  characterize  the  loss  be¬ 
havior.  Nevertheless,  the  formula  for  both  macrobending  and 
microbending  fiber  loss  employs  many  identical  mathematics. 


We  started  our  computer  model  effort  by  working  on  the  math¬ 
ematical  programming  of  the  constant  curvature  bending  loss  of 
step-index  single  mode  fiber.  Marcuse  has  shown  that  the  bending 
loss,  a  ,  can  be  expressed  in  terms  of  the  fiber  index  profile  and 
the  bend  radius  R  as:  (Ref. 47  in  Appendix  B) 
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a  =  finer  exp  [-  j  (^)  ft] 

2  73£  V2^(X-i  (7a)  K-.v73)] 


where 


«•  ( nc2  k2  -(J2)15  ,  7  =  (  (3 2-  nc£2  k2)^ 

V  =  «a("c2  -  V2)4  ,  K’Xa- 


a  is  the  fiber  cord  radius.  nc  and  are  the  refractive 

index  of  fiber  core  and  cladding  respectively.  X  is  optical  wave¬ 
length. 

A  computer  program  has  been  written  using  Professional  FOR¬ 
TRAN  as  its  source  language.  This  program  has  been  tested  on  an 
IBM  PC  AT  with  math  processor  and  should  run  on  IBM  PC,  XT,  or 
compatibles  with  a  math  coprocessor.  The  preliminary  program 
listing  is  included  as  Appendix  C.  The  program  calculates  the 
bending  loss  curves  as  a  function  of  fiber  parameters  and  bending 
radius  as  shown  in  Figure  3.  It  shows  that  the  bend  induced  loss 
depends  critically  on  the  fiber  core-cladding  refractive  index 
difference  and  the  bend  radius  R.  The  next  step  will  compare  the 
calculated  loss  values  with  measured  constant  curvature  bend  loss 
data  generated  from  fiber  samples  designed  for  use  in  high  speed 
missile  payout  dispensers.  Expected  discrepancies  between  the  two 
sets  of  loss  data  will  be  analyzed  for  improving  the  bending  loss 
model  as  well  as  for  bend  loss  measurement. 
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Although  this  bending  loss  program  is  written  specifically 
for  a  step-index  single  mode  fiber,  it  is  believed  that  it  will  be 
applicable  to  other  single  mode  fibers  with  different  refractive 
index  profiles.  This  extension  will  be  required  to  establish  an 
"equivalent  step-index  fiber  mode  field  size"  for  other  fibers  by 
matching  their  evanescent  field  tails  in  cladding  region  against 
an  ideal  step-index  fiber.  Theoretical  analysis  and  fiber  output 
spot  size  measurements  for  sample  fibers  will  be  conducted  to 
validate  this  bending  loss  analysis  concept. 

One  potential  application  for  the  constant  curvature  bending 
loss  study  is  to  evaluate  the  fiber  excess  loss  while  the  fiber  is 
subjected  to  a  constant  speed  payout.  The  payout  peel  point  cur¬ 
vature  is  suspected  to  be  a  major  loss  contributor  in  the  fiber 
payout  process.  A  mechanical  model  analysis  on  the  peel  point 
curvature  in  terms  of  fiber  parameters  and  payout  conditions  is 
\  _ng  improved  on  a  separate  project.  The  calculated  peel  point 
curvature  will  then  be  used  in  the  bending  loss  computer  program 
to  predict  the  fiber  loss  during  the  payout. 

Another  analysis  effort  currently  under  way  involves  the  col¬ 
lection  of  optical  loss  data  on  bobbin  wound  fibers  and  experimen¬ 
tal  data  on  different  loss  measurement  techniques.  Existing  fiber 
loss  data  indicate  that  winding  loss  and  low  temperature  excess 
loss  of  bobbin  wound  fibers  are  both  bending  loss  in  nature. 
Winding  geometry  and  the  material  thermal-mechanical  properties  of 
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fiber  buffer  layer  have  been  identified  as  the  prime  factors  in 
the  loss  analysis.  Additional  modeling  is  needed  to  formulate 
thermal  mechanical  effects  in  a  bobbin  wound  fiber  pack.  The 
stress  profile  of  the  fiber  pack  will  then  be  translated  into 
microbending  parameters  and  used  for  a  fiber  loss  prediction. 


FUTURE  PLAN 

The  immediate  plan  is  to  expand  the  bending  loss  computer 
program  to  include  periodic  bend  loss  analysis.  The  periodic  bend 
loss  program  will  then  further  be  expanded  to  cover  the 
microbending  loss  analysis  that  will  integrate  the  loss  contribu¬ 
tions  from  an  ensemble  of  microbend  perturbations  in  different 
spatial  frequencies.  This  effort,  along  with  bend  loss  measure¬ 
ment  on  sample  fibers,  will  be  the  primary  task  in  the  remainder 
of  the  Basic  Program.  Unless  the  option  of  the  proposed  Optional 
Program  is  exercised,  a  final  report  for  the  Basic  Program  will  be 
prepared  to  cover  the  finding  of  this  study.  If  the  option  is  ex¬ 
ercised,  the  study  effort  will  be  continued  as  shown  in  Figure  4 
in  the  Optional  Program.  The  results  will  be  presented  in  the 
form  of  a  progress  report  as  specified  by  the  contract. 

The  critical  task  of  the  Optional  Program  is  to  devise  prac¬ 
tical  fiber  loss  chacterization  schemes  that  will  yield  relevant 
input  data  useful  for  the  bending  loss  model.  Preliminary  loss 
measurements,  including  Optical  Time  Domain  Ref lectometer  ( OTDR) 
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the  bending  loss  model.  The  discrepancy  analysis  will  be  analyzed 
to  provide  new  leads  for  the  improvement  of  the  bending  loss 
model.  Similar  procedure  will  then  be  expanded  to  deal  with  both 
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c  This  program  calculates  the  normalize  frequency,  V,  with 
c  given  values:  lambda ( wave  1  eng th ) ,  delta,  A,  and  N(index 
c  of  refraction).  Futhermore,  it  calculates  BetaA,  KapaA. 
c  and  GamaA  of  the  transcendental  problem  of  the  Bessei  and 
c  the  Mod i f i ed - Be s s e 1  functions.  All  the  values  are  calculated 
c  in  dou b 1 e - p rec 1 s  1  on  . 


c  All  the  variables  are  implicitly  declared  real  except  I  and  M. 
c  Flfty(50)  elements  are  reserved  -for  the  arraies. 

IMPLICIT  REAL*8(A-H. J-L.N-Z) 

DIMENSION  RoA(SO),  R(50),  Alpha(50),  Argumt(50),  Expnt(50). 
c  DeA 1 phaA ( 50  )  ,  NuAl phaA ( 50 ) .  AlphaA(50),  AlphaL(50) 

c  There  is  an  input  file  called  ''IP”  that  this  program  reads  its  data 
c  directly  from. 

OPEN  (UNIT-8, FILE-'  IP1, STATUS- ' OLD  1  ) 

READ (8,10)  lambda,  delta,  A,  N.  alnc 
DO  5  1=1,19 
R  E AO ( 8 , 2  0  )  RoA(I) 

5  CONTINUE 

CLOSE  (UNIT-8) 

10  FORMAT ( 9X , F4.2,12X.F5.3.8X.F4.2,8X.F4.2,10X.F6.1) 

20  FORMAT ( 6X , FlO  .  1  ) 

Pi-3 . 14159265359 

c  This  is  where  the  NORMALIZE  FREQUENCY,  V,  is  calculated. 

V  =  ((2.*Pi)/'lambda)*A*N*(DSQRT(2.*delta)) 

PRINT* ,  ’ V  =  ■  ,  V 

PRINT*  ,  'Inc  - '  , alnc 
M-0 

c  This  where  the  calculation  of  upper  and  lower  bound  of  the 
c  transcendental  problem  is  calculated. 

UpBd-N*( (2*Pi)/ lambda )*A 

LwBd-N  * ( l -de 1 ta ) * ( ( 2*P1 )/ lambda ) *A 

BetaA-UpBd 

PxUpBd-UpBd 

FxLwBd-LwBd 

50  KapaA-DSQRT (  ( PxUpBd *  *  2 ) - ( Be t a A  *  *  2  )  ) 

GamaA-DSQRT ( ( BetaA* *2 ) - ( FxLwBd**2  )  ) 

CALL  Jo  (KapaA,  FJo) 

CALL  J 1  (KapaA.  PJl) 

CALL  Ko  (GamaA,  FKo ) 

CALL  K 1  (GamaA.  FK1  ) 

F J»KapaA* ( PJl/FJo) 

9  FK-GamaA* ( FK1 /FKo ) 

£r ror-FJ- FK 


This  accuracy  can  be  alter  to  agproxnately  i.xiOE-12 
IF  (ABS(Error)  .  LT  .  0.00000001)  GOTO  200 
IF  (Error  . GT .  0.0)  aSign*l 
IF  (Error  .LT.  0.0)  »Sign*2 
M-M*  1 

IF  (M  EQ.  1)  GOTO  100 
IF  (mFlag  .ME.  mSign)  THEN 
Up3d=NuUpBd 
LwBd-NuLwBd 
Be  taA-UpBd 
END  I  F 

100  »Flag*aSign 

CALL  Bound  ( U  pBd ,  LwBd.  BetaA.  NuUpBd.  NuLwBd.  delta,  alnc) 
IF  (M  .GT.  50000)  GOTO  1 
GOTO  50 

200  00  250  1*1.19 

Fixed  Lower  Bound  <  BetaA  <  Fixed  Upper  Bound 

Arguat (  I  )  *  (  2  .  /  3  .  )*(((( GaaaA ) *  *  3 ) / ( ( BetaA ) *«2 ) ) • ( RoA< I ) ) ) 
EGaaaA*(DSQRT(GaaaA)  )**3 
Expnt( I)*0£XP(-Arguat(  I)  ) 

NuAl phaA (  I)*(DSQRT(P1)*(  (KapaA)**2)*Expnt(  I  )  ) 

DeAlphaA! I ) -4*£GamaA*(V**2) *(DSQRT(RoA(I)))*(FKl**2) 

A 1 phaA ( I )*NuAlphaA( I )/DeAiphaA( I  ) 

R ( I ) *RoA ( I ) /A 
Alpha( I )*AlphaA(  I  )/A 
L  * ( P i *R (  I)  )  /  2 

Calculation  of  AIphaL 

AlphaL ( I ) -( 1-EXP( -Alpha! I ) *L ) ) /Alpha (  I  ) 

250  CONTINUE 

260  FORMAT!  1  laabda  - '  . F5 . 2 . 4X  .  '  del  ta  -  ',F5.3.4X.'A  *'. 

C  F 5 . 2 , 4 X ,  '  N  *  '  . F5 . 2 . 4X.  ' Inc  -  '  . F5 .  1  ) 

265  FORMAT (  '  '  ) 

Final  result  is  printed  out  in  the  output  file  called  " FORT9 
and  by  viewing  “ FORT9 "  will  display  all  final  result(s). 

WR I TE (9,280)  laabda.  delta.  A.  N.  alnc 
WRITE(9.265) 

3eta-BetaA/A 

270  FORMAT ( '  NORMALIZE  FREQUENCY  V  =  ’.F2D.18) 

272  FORMAT!’  Fixed  Upper  Bound  -  .F20.17) 

274  FORMAT!'  Fixed  Lower  Bound  *  ’.F20.17) 

276  FORMAT! 13X BetaA  -  ’.F20.17) 

278  FORMAT! '  KapaA  -  '. F20 . 18 , 5X GaaaA  -  ’.F20.18) 

279  FORMAT!'  Beta  -'.F20.17) 

WR I TE (  9 . 27 0  )  V 
WRITE(9.279)  Beta 
WRITE! 9 . 265 ) 

WRITE(9.272)  FxUpBd 

WR I TE ( 9 . 2 7 4  )  FxLwBd 
WRITE(9.276)  BetaA 
WRITE ( 9 , 265  ) 

WRITE(9.278)  KapaA,  GaaaA 
WRITE ( 9 . 265  ) 

280  FORMAT!  'R  (Spool  Radius):',  '  Alpha  (Energy  Loss):', 

c  '  AlphaL! EnergyLoss/Length )  :  ’) 

WRITE!  9 . 280.) 

290  FORMAT! 3X , FI l .2,6X.E21.15,5X.E21.15) 

DO  300  1-1,19 

WRITE(9,290)  R(I),  Alpha(I).  AlphaL(I) 

300  CONTINUE 

1  PRINT*,  'Nuaber  of  loop  is'.M 


c  -  - - - 

c  THIS  SUBROUTINE  CALCULATES  BetaA  WITH  GIVEN  BOUNDARY 

c  - 

c 

SUBROUTINE  Bound  (UpBd,  LwBd,  BetaA.  NuUpBd.  NuLwBd ,  delta,  alnc) 
IMPLICIT  R  E  A  L  *  8  !  A  -  Z  ) 

XI -UpBd 
X2  *  LwBd 

Del-(Xl*X2)/'aInc 

Be  taA=BetaA-De  1 

NuLwBd-Be  CaA 

NuUpBd-BetaA-Oel 

RETURN 

END 

c  Following  subroutines  are  the  Bessel  functions  of  zeroth 
c  and  of  first  order. 

c  - 

c  THIS  SUBROUTINE  CALCULATES  THE  Jo(X)  BESSEL  FUNCTIONS 

c  - 

c 

SUBROUTINE  Jo  (KapaA.  FJo  ) 

IMPLICIT  REAL*8(A-Z) 

X-KapaA 

IF  (X  .LE.  3.0)  GOTO  100 
T-3 . 0/X 

F  *0 . 79788456»T* ( -0 . 000Q0077*T* ( -0 . 00552740*-T* (  -0 . 000095  1  2  *  T  * 

C  (  0 . 00137237  *-T"  (  -  0  .  00072805  •-  T  *  0  .  00014476  )  )  )  )  ) 

THETA-X-0 . 78539816»T“ ( -  0 . 04166397»T*( -0 . 00003954*T*( 0 . 00262573-T* 
C  (  -0.00054125*T*{-0.00029333*-T*0. 00013558)  )  )  )  ) 

Q-l ■ 0/DSQRT(X) 

FJo-Q*F*DCOS(THETA) 

RETURN 

100  T-X*X/9.0 

FJo-1.0-T*(2. 2499997 - T* ( 1 . 2 858208-T«(0.3163866-T*(0.0444479-T* 

C  ( 0 . 0039444-T*0 . 0002100  ))))  ) 

RETURN 

END 

c  - 

c  THIS  SUBROUTINE  CALCULATES  THE  J1(X)  BESSEL  FUNCTIONS 

c  - 

c 

SUBROUTINE  J1  (KapaA,  FJ1) 

IMPLICIT  REA  L “ 8 ( A  - Z ) 

X-KapaA 

IF  (X  .LE.  3.0)  GOTO  100 
T  -  3 . 0  /  X 

FI -0 . 79788438*T* (  0 . 00000156-T* (  0 . 01659667  *T* ( 0 . 00017 105-T* 

C  (-0. 0024951 l*T*(0.00113653-T*000020033)  )  )  )  ) 

THETA- X-2 . 35619449-T*( 0.1249961 2-T*(0. 00005650-T*( *0.00637879* T* 

C  (0.00074348*T»(0.00079824-T*0. 00029166))))) 

Q-l . 0/DSQRT( X) 

FJ1-Q*F1*DC0S(THETA) 

RETURN 

100  T-X*X/9 . 0 

PJ1»X*{0.3-T*(0.56249983-T*(0.2:093573-T*(0.03954289-T* 

C  (0.00443319-T*(0.00031761-T*0. 00001109)))))) 

RETURN 

END 


Following  subroutines  are  the  aod 1 f i ed - Bes se  1  function  of  zeroth 
and  of  first  order. 

THIS  SUBROUTINE  CALCULATES  THE  Ko(X)  MOD  I  F  I  ED - BES S E L  FUNCTION 


SUBROUTINE  Ko  (GamaA.  FKo) 

IMPLICIT  R£AL*8(A-Z) 

X^GamaA 

IP  (X  .GT.  3.75)  GOTO  100 
T=X*X/ (3.75*3.75) 

F  I  o “  1  . 0  +  T  * ( 3 . 5156229»T*(  3 .  0899424*T* (  1  2067492~T*(0  2659732-T* 

C  ( 0 . 0360768-T*0  0045813))))) 

GOTO  200 
100  T-3.75/X 

IF  (X  . GT .  85.0)  X-85 . 0 

FIo-DEXP(  X  )  /  DSQRT(  X)  *(  0  .  39894228«-T*(  0  .  01328592*T*  (  0  .  00225319<-T* 

C  (  -0 . 00157565  >T* <  0 . 009  1628  1  *T* (  -0  . 02057706-T* ( 0  02635537  -T" 

C  (-0.01647633+T*0. 00392377)))))))) 

GOTO  200 

200  IF  (X  LT.  2.0)  GOTO  300 
T  »  2 . 0  /  X 

F Ko  =  DEXP( -  X ) /DSQRT( X ) * ( 1  .  2  5  3314 14-T* ( -0 . 07 83 235 8 »T* ( 0 . 02 1 89568  -T' 
C  (  -0 . 01062446<-T*  (  0  .  00587872*T*  (  -  0  .  00251540‘T*0  00053208  )  )  )  )  )  ) 
RETURN 

300  T*0.25*X*X 

IF  (X  . LT .  1  E-30  )  X«1 . E  30 

FKo»-DL0G( 0 . 5*X) *FIo-0 . 57721566»T*( 0 . 42278420»T*( 0 . 23069756-T* 

C  ( 0 . 03488590»T* ( 0 . 00262698*T* ( 0  00010750*T*0  000007  40  )  )  )  )  ) 

RETURN 
END 


THIS  SUBROUTINE  CALCULATES  THE  K1(X)  MOD  I F I  ED  -  BESSEL  FUNCTIONS 


SUBROUTINE  K1  (Ga»aA.  FK1 ) 

IMPLICIT  REAL*8(A-Z) 

X-GanaA 

IF  (X  .GT.  3.75)  GOTO  100 
T*X*X/ (3.73*3.75) 

FI1>X*(0.5-T*(0, 878905 94  <-T*(  0.51498869*T*(0.  15084934 *T» 

C  (0.026S8733*T*(0. 00 301532*T*0. 00032411)))))) 

GOTO  200 
100  T-3.75/X 

IF  (X  .GT.  83.0)  X-83.0 

FI 1*DEXP( X) /DSQRT( X) ■( 0 . 39894 2 2 8*T* ( -0 . 0398802 4 -T*( -0 . 0036201 8 -T* 
C  (0.00163801*T*(-0.01031535*T*(0.02282967*T*(-0.02895312-T* 

C  (0.01787654-T*0. 00420059)))))))) 

200  IP  (X  .LT.  2.0)  GOTO  300 
T-2 . 0/X 

FK1-0EXP( -X) /DSQRT( X) *( 1 . 2533 1 4 1 4 -T* ( 0 . 234986 19*T* ( -0 . 03655620 *T* 
C  (  0 . 01504268->-T*  (  -  0 . 00780353  +  T*  (  0  .  003256  14-T*0  .  00068245  )  )  )  )  )  ) 
RETURN 

300  T*0 . 25*X*X 

IF  (X  .LT.  l.E-30)  X  =  1  .E-30 

FK1-DL0G(  0 . 5*X)  *FI  1  *<  1  .  0*T*(  0  .  15443144<-T*(-0.67278579*T* 

C  (-0.18136897*T*(-0. 0191 9402 +  T«(-0.00110404-T*0. 00004686)))))) /X 
RETURN 
END 


■.  a  V.  a  V 


